# ************************************* #
# Author: Piatinskii M., Azov-black sea bench of VNIRO
# &copy 2019
# **************************************#
# make cleanup
rm(list = ls(all.names=TRUE))
# load libraries
library("FLCore")
library("FLXSA")
library("FLAssess")
library("ggplotFL")
# library("FLBRP") # no support for R 4.0, require compilation from *nix
library("FLash")
#install.packages("devtools")
#install_github("ices-tools-prod/msy")
library("msy")
library("rmarkdown")
library("dplyr")
library("magrittr")
library("tibble")
library("icesAdvice")

# ****** GLOBAL CONFIGURATIONS SECTION *******#

# stock object configurations
config.stock.name <- "Sprattus sprattus" # fish title
config.stock.desc <- "Sprat stock information in FLCore format from 2000years" # fish stock description
config.stock.fbar.min <- 1 # minimum fish age become fishing (Fbar min)
config.stock.fbar.max <- 3 # maximum fish age at fishing (Fbar max)
config.stock.no_harvestspawn <- TRUE # Set TRUE if there is no information about fishing before spawning
config.stock.no_mspawn <- TRUE # Set TRUE if no information about natural mortality before spawning
config.stock.no_discards <- TRUE # Set TRUE if no separated information about landing/discards exist (only catch)
config.stock.window.start <- 1994 # model start
config.stock.window.end <- 2022 # model end

# configure surveys indices 
config.index1.name <- "Sprat crimea cpue index" 
config.index1.desc <- "Fishery index information about sprat abundance at Crimea shelf from real fishing surveys"
config.index1.startf <- 0.3 # start month of survey in format month/12, example: 1 feb = 2/12 = 0.166
config.index1.endf <- 0.7 # end month of survey in format month/12
config.index1.type <- "number" # type of survey, allowed: number, biomass

# same with previous
config.index2.name <- "Sprat caucasian cpue index" 
config.index2.desc <- "Fishery index information about sprat abundance at Caucasian shelf from real fishing surveys"
config.index2.startf <- 0.3
config.index2.endf <- 0.7
config.index2.type <- "number"

# pre-analysis diagnostic procedures. Useful to choice the model parametrization, convergence, etc
config.iter.fse = TRUE # Fse choice procedure
config.iter.retro = TRUE # retrospective procedure
config.iter.fseage = TRUE # Fse at age impact procedure

# configure retrospective depth
config.retro.horizon <- 5

# Shrinkage standard error max value. Default value Fse = 0.5
# Higher shrinkage se leads to high unrestraint of estimating procedure. 
config.tuning.xsa <- FLXSA.control(
  maxit = 100, # maximum covergence iterations
  fse = 1.5, # shrinkage s.e value
  shk.n = T,
  shk.f = T,
  shk.yrs = 4,
  tsrange = 15, # weighting time series length (from latest year)
  tspower = 3, # weighting procedure sequence (1 = linear, 2 = square, 3 = cubic)
  qage = 2, # age after catchability not depends on age
)

# sr model configuration
config.srmodel.eqsim.use = TRUE # should we use eqsim package insted of FLSR default?
config.srmodel.use_geomean = FALSE # set FALSE to use EQSIM SR model results to forecasting

# eqsim model configurations (if config.srmodel.eqsim.use = TRUE)
config.srmodel.eqsim.models = c("Segreg") # types of model to fit and compare
config.srmodel.eqsim.simulations = 10 # number of simulations
config.srmodel.eqsim.shift_years <- 1994:2005 # remove years from sr-model fitting

# flsr model configuration (if config.srmodel.eqsim.use = FALSE)
config.srmodel.type <- segregAR1() # SR model type: ricker(), bevholt(), geomean(), segreg(), shepherd(), cushing(), AR in model name = AutoCorrelation fix
config.srmodel.optimization = FALSE # Use optimization procedure for SR model fitting?

# configure MSY reference point calculation: Blim, Bpa, Fmsy
config.brp.blim_multiplier <- 0.3 # multiplier of SSBmax to get Blim ref.point, allowed: 0.2 - 0.4
config.brp.trim_uncertainity <- TRUE # trim uncerteinity extreme deviations? Optimization procedure for Bpa estimate
config.brp.current_range <- 2005:2022 # shift previous population states from BRP calculations

# Custom forecast scenarious. By default performed: F_SQ, F_01, F_msy scenarious.
# There you can define f(y+1), f(y+2), f(y+3) targets for 2 scenarious
# precautionary scenario: Fpa = (B-Blim/Bmsy-Blim) * Fmsy = 0.51
config.forecast.scenario1.f_year1 <- 0.51
config.forecast.scenario1.f_year2 <- 0.51
config.forecast.scenario1.f_year3 <- 0.51
config.forecast.scenario1.f_year4 <- 0.51
config.forecast.scenario1.f_year5 <- 0.51

# Scenario MSY at M = 0,95
# by old STECF (2010) results
config.forecast.scenario2.f_year1 <- 0.91
config.forecast.scenario2.f_year2 <- 0.91
config.forecast.scenario2.f_year3 <- 0.91
config.forecast.scenario2.f_year4 <- 0.91
config.forecast.scenario2.f_year5 <- 0.91

# additional scenario for manual Fmsy and F_01 estimates
config.forecast.manual_F01 <- 0.605
config.forecast.manual_FMSY <- 0.64

# custom calculations. Calculate partial recruitment for Y/R analysis
config.custom.partial.years <- 2004:2020

# custom calculations. Calculate F01 reference point from Y/R analysis
config.custom.f01.years <- 2004:2020 # years used to calculate mean(wt), mean(Fi/Fmax), mean(M)
config.custom.f01.plusgroup = TRUE # is last group is plus group?
config.custom.f01.oldest = 5 # oldest fish found in waters (not in stock, that means "even somewhen exist in sea")
# ******** GLOBAL CONFIGURATION DONE!!! **** #

# ********************************************************* #
# Do not change enything below if you have no idea what is going on
# ********************************************************* #

# simple write functions
savePic <- function(path, imgObj, width=1900, height=1600, res=200) {
  png(path, width=width, height=height, res=res, type="cairo")
  print(plot(imgObj) + theme_bw())
  dev.off()
}
saveText <- function(path, object) {
  con <- file(path, encoding = "UTF-8")
  sink(con)
  print(object)
  sink()
  close(con)
}

# get mbar from stock object
mbar <- function(stock) {
  stock.min <- as.numeric(stock@range["min"])
  stock.max <- as.numeric(stock@range["max"])
  
  bar.min <- as.numeric(stock@range["minfbar"])
  bar.max <- as.numeric(stock@range["maxfbar"])
  
  diff.min <- 1 + bar.min - stock.min
  diff.max <- 1 + bar.max - bar.min
  if (diff.max > stock.max+1)
    diff.max <- stock.max+1
  
  vec <- (m(stock)[diff.min:diff.max,])
  bar <- apply(vec, 2, mean)
  return(bar)
}



# read input data - stock
tmp.stock <- read.csv("input/stock.csv", sep=",")
data.stock <- as.FLStock(tmp.stock)
# read input data - index
tmp.idx1.idx <- read.csv("input/survey1/index.csv", sep=",")
tmp.idx1.eff <- read.csv("input/survey1/effort.csv", sep=",")
tmp.idx2.idx <- read.csv("input/survey2/index.csv", sep=",")
tmp.idx2.eff <- read.csv("input/survey2/effort.csv", sep=",")

data.index1 <- FLIndex(index=as.FLQuant(tmp.idx1.idx), effort=as.FLQuant(tmp.idx1.eff), catch.n=as.FLQuant(tmp.idx1.idx))
data.index2 <- FLIndex(index=as.FLQuant(tmp.idx2.idx), effort=as.FLQuant(tmp.idx2.eff), catch.n=as.FLQuant(tmp.idx2.idx))

remove(tmp.stock, tmp.idx1.idx, tmp.idx1.eff, tmp.idx2.idx, tmp.idx2.eff)
print("[XSA toolkit]: Load data: OK")
# Set name,desc from cfgs
name(data.stock) <- config.stock.name
desc(data.stock) <- config.stock.desc
# set fbar ranges
range(data.stock)[c("minfbar", "maxfbar")] <- c(config.stock.fbar.min, config.stock.fbar.max)

# shift data starting from 2000
data.stock <- window(data.stock, start = config.stock.window.start, end = config.stock.window.end)

# set harvest.spwn and m.spwn = 0 if no data exist (or XSA throw error)
if (config.stock.no_harvestspawn) {
  harvest.spwn(data.stock) <- 0
}
if (config.stock.no_mspawn) {
  m.spwn(data.stock) <- 0
}
# if no info about discards/landings separatly than set discards = 0, landings = catch
if (config.stock.no_discards) {
  discards(data.stock) <- 0
  discards.n(data.stock) <- 0
  discards.wt(data.stock) <- catch.wt(data.stock)
  landings(data.stock) <- catch(data.stock)
  landings.n(data.stock) <- catch.n(data.stock)
  landings.wt(data.stock) <- catch.wt(data.stock)
}
stock.wt(data.stock) <- catch.wt(data.stock)

# set cfgs for index1 and index2
type(data.index1) <- config.index1.type
name(data.index1) <- config.index1.name
desc(data.index1) <- config.index1.desc
range(data.index1)[c("startf", "endf")] <- c(config.index1.startf, config.index1.endf)

type(data.index2) <- config.index2.type
name(data.index2) <- config.index2.name
desc(data.index2) <- config.index2.desc
range(data.index2)[c("startf", "endf")] <- c(config.index2.startf, config.index2.endf)
print("[XSA toolkit]: FLStock and FLIndex object configured successful")

data.indices <- FLIndices(data.index1, data.index2)

# make Fse comparision test
if (config.iter.fse) {
  print("[XSA toolkit]: Fse impact factor estimate")
  steps <- seq(0.5, 2.5, by = 0.5)
  fseIters <- c()
  for (i in 1:length(steps)) {
    ctl <- config.tuning.xsa
    ctl@fse <- steps[i]
    res <- FLXSA(data.stock, data.indices, ctl) + data.stock
    fseIters <- c(fseIters, res)
  }
  
  result.fse.summary <- FLStocks(fseIters)
  names(result.fse.summary) <- steps
  # save output to pic
  savePic("output/1_FseChoose/summary_all.png", result.fse.summary)
  # get latest 5 year window
  lastYear <- result.fse.summary[[1]]@range["maxyear"]
  tmpObj <- window(result.fse.summary, lastYear-config.retro.horizon, lastYear)
  savePic("output/1_FseChoose/summary_short.png", tmpObj)
  remove(lastYear, tmpObj)
  # save text output - ssb
  con <- file("output/1_FseChoose/ssb.txt", encoding = "UTF-8")
  sink(con)
  
  for (i in 1:length(steps)) {
    print(sprintf("============= SSB at Fse = %s", steps[i]))
    print(ssb(result.fse.summary[[i]]))
    cat("\r\n")
  }
  
  sink()
  close(con)
  # save text output - rec
  con <- file("output/1_FseChoose/rec.txt", encoding = "UTF-8")
  sink(con)
  
  for (i in 1:length(steps)) {
    print(sprintf("============= Recruitment at Fse = %s", steps[i]))
    print(rec(result.fse.summary[[i]]))
    cat("\r\n")
  }
  
  sink()
  close(con)
  # save text output - fbar
  con <- file("output/1_FseChoose/fbar.txt", encoding = "UTF-8")
  sink(con)
  
  for (i in 1:length(steps)) {
    print(sprintf("============= Fbar at Fse = %s", steps[i]))
    print(fbar(result.fse.summary[[i]]))
    cat("\r\n")
  }
  
  sink()
  close(con)
  
  remove(steps, fseIters, ctl, res, con)
}

if (config.iter.retro) {
  print("[XSA toolkit]: Retrospective analysis")
  # make retrospective analysis
  result.retro.summary <- list()
  endYear <- data.stock@range["maxyear"]
  startYear <- endYear-config.retro.horizon
  years <- startYear:endYear
  frange <- seq(0.1, 2.5, by=0.1)
  for (i in frange) {
    res <- tapply(years, 1:length(years), function(x) {
      ctl <- config.tuning.xsa
      ctl@fse <- i
      window (data.stock, end=x) + FLXSA(window(data.stock, end=x), data.indices, ctl)
    })
    res <- FLStocks(res)
    res@names <- ac(c(years))
    result.retro.summary[[ac(i)]] <- res
    
    # save output as img
    savePic(sprintf("output/2_RetroFse/%s.png", i), res)
    remove(res)
  }
}

print("[XSA toolkit]: Retrospective Mohn's rho tests")
# compare XSA estimates in terminal year by -5 year retrospective
lastYear <- as.numeric(range(data.stock)["maxyear"])
ssb.base <- as.vector(ssb(window(data.stock+FLXSA(data.stock, data.indices, config.tuning.xsa), start=lastYear-5)))
f.base <- as.vector(fbar(window(data.stock+FLXSA(data.stock, data.indices, config.tuning.xsa), start=lastYear-5)))
ssb.retro <- list()
f.retro <- list()
for (i in 1:config.retro.horizon) { # prepare retrospective matrix for mohn's rho test
  # make retro
  tmp.xsa <- FLXSA(window(data.stock, end=(lastYear-i)), window(data.indices, end=(lastYear-i)), config.tuning.xsa)
  tmp.xsa.st <- data.stock + tmp.xsa
  
  ssb.t <- as.vector(ssb(window(tmp.xsa.st, start=lastYear-config.retro.horizon)))
  f.t <- as.vector(fbar(window(tmp.xsa.st, start=lastYear-config.retro.horizon)))
  
  
  ssb.retro[[i]] <- ssb.t
  f.retro[[i]] <- f.t
}

result.retro.diagnostic.ssb <- data.frame(ssb.base, ssb.retro[[1]], ssb.retro[[2]], ssb.retro[[3]], ssb.retro[[4]], ssb.retro[[5]])
rownames(result.retro.diagnostic.ssb) <- (lastYear-config.retro.horizon):lastYear
result.retro.diagnostic.f <- data.frame(f.base, f.retro[[1]], f.retro[[2]], f.retro[[3]], f.retro[[4]], f.retro[[5]])
rownames(result.retro.diagnostic.f) <- (lastYear-config.retro.horizon):lastYear
names(result.retro.diagnostic.ssb) <- c("ssb.base", "-1", "-2", "-3", "-4", "-5")
names(result.retro.diagnostic.f) <- c("f.base", "-1", "-2", "-3", "-4", "-5")

# replace NaN's (x/0) to NA (not available)
is.nan.data.frame <- function(x) {
  return (do.call(cbind, lapply(x, is.nan)))
}

result.retro.diagnostic.f[is.nan(result.retro.diagnostic.f)] <- NA
result.retro.diagnostic.ssb[is.nan(result.retro.diagnostic.ssb)] <- NA

result.retro.mohnrho.ssb <- mohn(result.retro.diagnostic.ssb, peels = config.retro.horizon)
result.retro.mohnrho.fbar <- mohn(result.retro.diagnostic.f, peels = config.retro.horizon)
remove(ssb.base, f.base, ssb.retro, f.retro)

print(paste0("Rho ssb=", round(result.retro.mohnrho.ssb, 3)))
print(paste0("Rho F=", round(result.retro.mohnrho.fbar, 3)))

# calculate Fse impact to different age groups
if (config.iter.fseage) {
  print("[XSA toolkit]: Fse impact to F-at-age estimates")
  steps <- seq(0.5, 2.5, by = 0.5)
  f.res <- propagate(harvest(data.stock), length(steps))
  for (i in 1:length(steps)) {
    raw.control <- config.tuning.xsa
    raw.control@fse <- steps[i]
    iter(f.res, i) <- harvest(FLXSA(data.stock, data.indices, raw.control))
  }
  endYear <- data.stock@range["maxyear"]
  startYear <- endYear - 9 
  tmpObj <- xyplot(data ~ year | age, groups = iter, data = f.res, ylab="F", type = "l", xlim = c(startYear:endYear))
  savePic("output/3_FseImpactFAtAge/summary_all.png", tmpObj)
  remove(startYear, endYear, steps, f.res, raw.control, tmpObj)
}

# save survey regression age-vs-age test 
suppressWarnings(savePic("output/4_SurveysRegressions/summary_all.png", data.indices))

suppressWarnings(savePic("output/4_SurveysRegressions/survey_1.png", data.indices[[1]]))
suppressWarnings(savePic("output/4_SurveysRegressions/survey_2.png", data.indices[[2]]))

print("[XSA toolkit]: Performing XSA analysis...")
# fit xsa
result.xsa.fit <- FLXSA(data.stock, data.indices, config.tuning.xsa)
# build final stock + xsa estimates
result.xsa.stock <- data.stock + result.xsa.fit

# estimate uncertainty right now!
# get f vals from xsa results
tmp.fvals <- as.vector(fbar(result.xsa.stock))
# get whiskers statistics for distribution
tmp.wstats <- boxplot.stats(tmp.fvals)$stats
# find outliers and remove from fvals
tmp.fvals <- tmp.fvals[tmp.fvals > tmp.wstats[[1]]]
tmp.fvals <- tmp.fvals[tmp.fvals < tmp.wstats[[5]]]
# get fixed uncerteinity estimation
result.xsa.sigma <- sd(tmp.fvals)
remove(tmp.fvals, tmp.wstats)

# save results 
endYear <- result.xsa.stock@range["maxyear"]
F_SQ = mean(fbar(window(result.xsa.stock, (endYear-2), endYear)))
savePic("output/5_XsaEstimates/summary_all.png", result.xsa.stock)
savePic("output/5_XsaEstimates/stock_numbers_at_age.png", stock.n(result.xsa.stock))
savePic("output/5_XsaEstimates/f_at_age.png", harvest(result.xsa.stock))

saveText("output/5_XsaEstimates/ssb.txt", ssb(result.xsa.stock))
saveText("output/5_XsaEstimates/rec.txt", rec(result.xsa.stock))
saveText("output/5_XsaEstimates/fbar.txt", fbar(result.xsa.stock))
saveText("output/5_XsaEstimates/f.txt", harvest(result.xsa.stock))
saveText("output/5_XsaEstimates/stock.n.txt", stock.n(result.xsa.stock))
saveText("output/5_XsaEstimates/diagnostics.txt", diagnostics(result.xsa.fit))

# investigate uncertainty factor (build conf.intervals for plots)
pics <- FLQuants(
  Rec = rnorm(1000, rec(result.xsa.stock), result.xsa.sigma * rec(result.xsa.stock)),
  SSB = rnorm(1000, ssb(result.xsa.stock), result.xsa.sigma * ssb(result.xsa.stock)),
  F = rnorm(1000, fbar(result.xsa.stock), result.xsa.sigma * fbar(result.xsa.stock)),
  Catch = catch(result.xsa.stock)
)

savePic("output/5_XsaEstimates/summary_uncertainty.png", pics)
remove(endYear, pics)

# perform residual diagnostics
pic <- bubbles(age~year, data=result.xsa.fit@index.res, warning=FALSE)
suppressWarnings(savePic("output/6_Residuals/summary_all.png", pic))
remove(pic)

pic <- bubbles(age~year, data=result.xsa.fit@index.res[1], warning=FALSE)
suppressWarnings(savePic("output/6_Residuals/summary_survey1.png", pic))
remove(pic)

pic <- bubbles(age~year, data=result.xsa.fit@index.res[2], warning=FALSE)
suppressWarnings(savePic("output/6_Residuals/summary_survey2.png", pic))
remove(pic)

pic <- xyplot(
  data ~ year | ac(age) + qname, 
  data=index.res(result.xsa.fit), 
  panel = function(x,y, ...){
    panel.xyplot(x, y, ...)
    panel.loess(x,y, ...)
    panel.abline(h=0, col="gray", lty=2)
  }, 
  main="Surveys log catchability residuals at age")
suppressWarnings(savePic("output/6_Residuals/catch_at_age_summary.png", pic))
remove(pic)

con <- file("output/6_Residuals/residuals.txt", encoding = "UTF-8")
sink(con)
print(result.xsa.fit@index.res)
cat("\r\n")
print("===== Summary statistics")
print(summary(result.xsa.fit@index.res))
sink()
close(con)

# analyse recruits-per-ssb unit 
result.sr.recssb <- rec(result.xsa.stock) / ssb(result.xsa.stock)
tmp <- as.data.frame(log(result.sr.recssb))
png("output/7_SR_BRP/rec_per_ssb.png", width=1900, height=1600, res=200, type="cairo")
print(plot(tmp$year, tmp$data, type="b", xlab="Year", ylab = "Ln(Recruit/SSB)"))
dev.off()
remove(tmp)


# calculate stock-recruitment relation
# be careful to use bevholt or ricker approximation: there is different approach between eqsim and flsr packages!!!
# Bevholt (eqsim): rec ~ log(a * ssb / (1 + b*ssb))
# Bevholt (flsr): rec ~ (a * ssb / (b+ssb))
# ricker (eqsim): log(a) + log(b) - b*ssb
# ricker (flsr): a * ssb * exp(-b * ssb)
# segreg(flsr,eqsim): ifelse(ssb < b, ssb*b, a*b)
if (config.srmodel.eqsim.use) {
  # build SR model based on EQSIM package and MSY approach
  print("[XSA toolkit]: Built SR model by EQSIM")
  
  result.eqsim.srmodel <- eqsr_fit(result.xsa.stock, 
                                   nsamp = config.srmodel.eqsim.simulations, 
                                   models = config.srmodel.eqsim.models,
                                   remove.years = config.srmodel.eqsim.shift_years
  )
  
  savePic("output/7_SR_BRP/eqsim_model.png", eqsr_plot(result.eqsim.srmodel, ggPlot = TRUE))
  saveText("output/7_SR_BRP/eqsim_model.txt", result.eqsim.srmodel$sr.det)
} else {
  # built Stock-Recruitment model based on FLSR and fmle()
  print("[XSA toolkit]: Built SR model by fmle()")
  result.sr.stock <- as.FLSR(window(result.xsa.stock, 
                                    start = config.stock.window.start, 
                                    end = config.stock.window.end))
  model(result.sr.stock) <- config.srmodel.type
  if (config.srmodel.optimization) {
    print("[XSA toolkit]: Performing SANN (Nelder-Mead, quasi-Newton, etc) optimization SR model. This can take a lot of time ...")
    result.sr.model <- fmle(result.sr.stock, method="L-BFGS-B", control=list(trace=0))
  } else {
    result.sr.model <- fmle(result.sr.stock, control=list(trace=0))
  }
  # save output
  png("output/7_SR_BRP/flsr_model.png", width=1900, height=1600, res=200, type="cairo")
  suppressWarnings(suppressMessages(plot(result.sr.model)))
  dev.off()
  
  png("output/7_SR_BRP/flsr_profile.png", width=1900, height=1600, res=200, type="cairo")
  suppressWarnings(profile(result.sr.model, plot=TRUE, trace = FALSE))
  dev.off()
  
  saveText("output/7_SR_BRP/flsr_model.txt", summary(result.sr.model)) 
}

print("[XSA toolkit]: Calculate BRP")

# calculate Bvirgin in assumtion of SSBmax + Ci
#startYear <- as.numeric(data.stock@range["minyear"])
#terminalYear <- as.numeric(data.stock@range["maxyear"])
startYear <- min(config.brp.current_range)
terminalYear <- max(config.brp.current_range)
df <- data.frame(year = startYear:terminalYear, ssb = as.vector(ssb(window(result.xsa.stock, start = startYear, end = terminalYear))), catch = as.vector(catch(window(result.xsa.stock, start = startYear, end = terminalYear))))
tmp.maxssb <- df[df$ssb == max(df$ssb),]
tmp.biomass <- tmp.maxssb$ssb + tmp.maxssb$catch

# process manual BRP points
result.brp.blim <- tmp.biomass * config.brp.blim_multiplier
# get Bpa based on Blim and uncerteinity sd
result.brp.bpa <- result.brp.blim * exp(qnorm(0.95) * result.xsa.sigma)

# estimate BRP from EQSIM
lastYear <- as.numeric(result.xsa.stock@range["maxyear"])
if (config.srmodel.eqsim.use) {
  result.eqsim.brp <- eqsim_run(result.eqsim.srmodel, 
                               bio.years = c(lastYear-4, lastYear),
                               bio.const = TRUE,
                               extreme.trim = c(0.1, 0.9), 
                               Bpa = result.brp.bpa,
                               Blim = result.brp.blim) 
  result.brpsr.blim <- as.numeric(result.eqsim.srmodel$sr.det["b"])
  result.brpsr.bpa <- as.numeric(result.eqsim.srmodel$sr.det["b"]*exp(qnorm(0.95)*result.xsa.sigma))
}

if (config.srmodel.use_geomean == TRUE) {
  # create geomean SR model by latest 5 years geomean recruitment numbers
  tmp.rec.mean <- exp(mean(log(window(rec(result.xsa.stock), lastYear-5, lastYear-1))))
  # build geomean SR model
  result.sr.model <- FLSR(model = "geomean", params=FLPar(tmp.rec.mean))
  remove(tmp.rec.mean, tmp.maxssb, df, tmp.biomass)
} else {
  # todo: process model
  result.sr.model <- FLSR(model = "segreg")
  result.sr.model@params["a"] <- as.numeric(result.eqsim.srmodel$sr.det["a"])
  result.sr.model@params["b"] <- as.numeric(result.eqsim.srmodel$sr.det["b"])
}

# plot SSB graphic with reference points Blim, Bpa
minYear <- as.numeric(result.xsa.stock@range["minyear"])
png("output/7_SR_BRP/ssb_blim_bpa_manual.png", width=1900, height=1600, res=200, type="cairo")
print(
plot(ssb(result.xsa.stock)) + 
  geom_hline(aes(yintercept=result.brp.blim), linetype=2, col="red") + 
  annotate(geom = "text", x = minYear+1, y = result.brp.blim+result.brp.blim*0.05, label = sprintf("Blim=%s", round(result.brp.blim, 2)), color = "red") +
  geom_hline(aes(yintercept=result.brp.bpa), linetype=2, col="blue") +
  annotate(geom = "text", x = minYear+1, y = result.brp.bpa + result.brp.bpa*0.05, label=sprintf("Bpa=%s", round(result.brp.bpa, 2)), color="blue") + 
  theme_bw()
)
      
dev.off()
remove(minYear)

# plot SSB with reference points by SR-model
# @todo
minYear <- as.numeric(result.xsa.stock@range["minyear"])
png("output/7_SR_BRP/ssb_blim_bpa_srmodel.png", width=800, height=600, res=150, type="cairo")
print(
  plot(ssb(result.xsa.stock)) + 
    geom_hline(aes(yintercept=result.brpsr.blim), linetype=2, col="red") + 
    annotate(geom = "text", x = minYear+2, y = result.brpsr.blim*1.2, label = sprintf("Blim=%s", round(result.brpsr.blim, 2)), color = "red") +
    geom_hline(aes(yintercept=result.brpsr.bpa), linetype=2, col="blue") +
    annotate(geom = "text", x = minYear+2, y =result.brpsr.bpa*1.2, label=sprintf("Bpa=%s", round(result.brpsr.bpa, 2)), color="blue") + 
    theme_bw()
)

dev.off()
remove(minYear)

# plot Fbar with F0.1 reference point
minYear <- as.numeric(result.xsa.stock@range["minyear"])
png("output/7_SR_BRP/fbar_f01_manual.png", width=1900, height=1600, res=200, type="cairo")
print(plot(fbar(result.xsa.stock)) + 
        geom_hline(aes(yintercept=config.forecast.manual_F01), linetype=2, col="red") + 
        annotate(geom="text", x = minYear, y = config.forecast.manual_F01+0.1, label = sprintf("F0.1=%s", round(config.forecast.manual_F01,3)), color="red") + 
        theme_bw()
      )

dev.off()
remove(minYear)

# plot E = 0.4 BRP reference point (ICES advice)
# E = F/(F+M)
# E = exploitation rate (*100 = percentage of usage)
result.er.table <- fbar(result.xsa.stock) / (fbar(result.xsa.stock) + mbar(result.xsa.stock))
png("output/7_SR_BRP/e_rate04.png", width=1900, height=1600, res=200, type="cairo")
print(plot(result.er.table) + ylab(label = "E = F/Z") + 
  geom_hline(yintercept = 0.4, color = "red", linetype="dashed") + 
  annotate(geom = "text", x = as.numeric(result.xsa.stock@range["minyear"]), y = 0.41, label="E = 0.4", color = "red") + 
  theme_bw())
dev.off()

# Prepare forecast procedure
print("[XSA toolkit]: Forecast procedure with different scenarious")
data.stf <- stf(result.xsa.stock, nyears=5)
lastYear <- as.numeric(result.xsa.stock@range["maxyear"])

# perform 3-year short-term forecasts by F_SQ
print("[XSA toolkit]: Short-term forecast by F_SQ")
target.tmp.F_SQ <- data.frame(year=(lastYear+1):(lastYear+5), quantity="f", val=F_SQ)
target.forecast.F_SQ <- fwdControl(target.tmp.F_SQ)
result.forecast.F_SQ <- fwd(data.stf, sr=result.sr.model, control=target.forecast.F_SQ)
savePic("output/8_Forecasts/F_SQ/summary_all.png", result.forecast.F_SQ)
savePic("output/8_Forecasts/F_SQ/summary_short.png", window(result.forecast.F_SQ, lastYear, lastYear+5))

saveText("output/8_Forecasts/F_SQ/ssb.txt", ssb(result.forecast.F_SQ))
saveText("output/8_Forecasts/F_SQ/rec.txt", rec(result.forecast.F_SQ))
saveText("output/8_Forecasts/F_SQ/fbar.txt", fbar(result.forecast.F_SQ))
saveText("output/8_Forecasts/F_SQ/f.txt", harvest(result.forecast.F_SQ))
saveText("output/8_Forecasts/F_SQ/stock.n.txt", stock.n(result.forecast.F_SQ))
saveText("output/8_Forecasts/F_SQ/ssb.txt", ssb(result.forecast.F_SQ))
saveText("output/8_Forecasts/F_SQ/catch.txt", catch(result.forecast.F_SQ))

remove(target.tmp.F_SQ)

# perform 3-year short-term forecasts by manual scenarius F_msy
print("[XSA toolkit]: Manual short-term forecast at F_MSY")
target.tmp.F_01manual <- data.frame(year=(lastYear+1):(lastYear+5), quantity="f", val=config.forecast.manual_F01)
target.forecast.F_01manual <- fwdControl(target.tmp.F_01manual)
result.forecast.F_01manual <- fwd(data.stf, sr=result.sr.model, control=target.forecast.F_01manual)

savePic("output/8_Forecasts/F_01manual/summary_all.png", result.forecast.F_01manual)
savePic("output/8_Forecasts/F_01manual/summary_short.png", window(result.forecast.F_01manual, lastYear, lastYear+5))

saveText("output/8_Forecasts/F_01manual/ssb.txt", ssb(result.forecast.F_01manual))
saveText("output/8_Forecasts/F_01manual/rec.txt", rec(result.forecast.F_01manual))
saveText("output/8_Forecasts/F_01manual/fbar.txt", fbar(result.forecast.F_01manual))
saveText("output/8_Forecasts/F_01manual/f.txt", harvest(result.forecast.F_01manual))
saveText("output/8_Forecasts/F_01manual/stock.n.txt", stock.n(result.forecast.F_01manual))
saveText("output/8_Forecasts/F_01manual/ssb.txt", ssb(result.forecast.F_01manual))
saveText("output/8_Forecasts/F_01manual/catch.txt", catch(result.forecast.F_01manual))

# perform 3-year short-term forecasts by manual scenarius F_MSY=config
print("[XSA toolkit]: Manual short-term forecast at F_MSY")
target.tmp.F_MSY <- data.frame(year=(lastYear+1):(lastYear+5), quantity="f", val=config.forecast.manual_FMSY)
target.forecast.F_MSY <- fwdControl(target.tmp.F_MSY)
result.forecast.F_MSY <- fwd(data.stf, sr=result.sr.model, control=target.forecast.F_MSY)

savePic("output/8_Forecasts/F_MSY/summary_all.png", result.forecast.F_MSY)
savePic("output/8_Forecasts/F_MSY/summary_short.png", window(result.forecast.F_MSY, lastYear, lastYear+5))

saveText("output/8_Forecasts/F_MSY/ssb.txt", ssb(result.forecast.F_MSY))
saveText("output/8_Forecasts/F_MSY/rec.txt", rec(result.forecast.F_MSY))
saveText("output/8_Forecasts/F_MSY/fbar.txt", fbar(result.forecast.F_MSY))
saveText("output/8_Forecasts/F_MSY/f.txt", harvest(result.forecast.F_MSY))
saveText("output/8_Forecasts/F_MSY/stock.n.txt", stock.n(result.forecast.F_MSY))
saveText("output/8_Forecasts/F_MSY/ssb.txt", ssb(result.forecast.F_MSY))
saveText("output/8_Forecasts/F_MSY/catch.txt", catch(result.forecast.F_MSY))


# short-term forecast at scenario1
print("[XSA toolkit]: Short-term forecast by Scenario1")
target.tmp.F_scen1 <- data.frame(year=(lastYear+1):(lastYear+5), quantity="f", val=c(
  config.forecast.scenario1.f_year1,
  config.forecast.scenario1.f_year2,
  config.forecast.scenario1.f_year3,
  config.forecast.scenario1.f_year4,
  config.forecast.scenario1.f_year5
))
target.forecast.Fsc1 <- fwdControl(target.tmp.F_scen1)
result.forecast.F_scenario1 <- fwd(data.stf, sr=result.sr.model, control=target.forecast.Fsc1)

savePic("output/8_Forecasts/F_scenario1/summary_all.png", result.forecast.F_scenario1)
savePic("output/8_Forecasts/F_scenario1/summary_short.png", window(result.forecast.F_scenario1, lastYear, lastYear+5))

saveText("output/8_Forecasts/F_scenario1/ssb.txt", ssb(result.forecast.F_scenario1))
saveText("output/8_Forecasts/F_scenario1/rec.txt", rec(result.forecast.F_scenario1))
saveText("output/8_Forecasts/F_scenario1/fbar.txt", fbar(result.forecast.F_scenario1))
saveText("output/8_Forecasts/F_scenario1/f.txt", harvest(result.forecast.F_scenario1))
saveText("output/8_Forecasts/F_scenario1/stock.n.txt", stock.n(result.forecast.F_scenario1))
saveText("output/8_Forecasts/F_scenario1/ssb.txt", ssb(result.forecast.F_scenario1))
saveText("output/8_Forecasts/F_scenario1/catch.txt", catch(result.forecast.F_scenario1))

remove(target.tmp.F_scen1)

# st forecast scenario2
print("[XSA toolkit]: Short-term forecast by Scenario2")
target.tmp.F_scen2 <- data.frame(year=(lastYear+1):(lastYear+5), quantity="f", val=c(
  config.forecast.scenario2.f_year1,
  config.forecast.scenario2.f_year2,
  config.forecast.scenario2.f_year3,
  config.forecast.scenario2.f_year4,
  config.forecast.scenario2.f_year5
))
target.forecast.Fsc2 <- fwdControl(target.tmp.F_scen2)
result.forecast.F_scenario2 <- fwd(data.stf, sr=result.sr.model, control=target.forecast.Fsc2)

savePic("output/8_Forecasts/F_scenario2/summary_all.png", result.forecast.F_scenario2)
savePic("output/8_Forecasts/F_scenario2/summary_short.png", window(result.forecast.F_scenario2, lastYear, lastYear+5))

saveText("output/8_Forecasts/F_scenario2/ssb.txt", ssb(result.forecast.F_scenario2))
saveText("output/8_Forecasts/F_scenario2/rec.txt", rec(result.forecast.F_scenario2))
saveText("output/8_Forecasts/F_scenario2/fbar.txt", fbar(result.forecast.F_scenario2))
saveText("output/8_Forecasts/F_scenario2/f.txt", harvest(result.forecast.F_scenario2))
saveText("output/8_Forecasts/F_scenario2/stock.n.txt", stock.n(result.forecast.F_scenario2))
saveText("output/8_Forecasts/F_scenario2/ssb.txt", ssb(result.forecast.F_scenario2))
saveText("output/8_Forecasts/F_scenario2/catch.txt", catch(result.forecast.F_scenario2))

remove(target.tmp.F_scen2)

# define stock names to summary plot
name(result.forecast.F_01manual) <- sprintf("F(0.1)=%s", config.forecast.manual_F01)
name(result.forecast.F_MSY) <- sprintf("F(MSY)=%s", config.forecast.manual_FMSY)
name(result.forecast.F_scenario1) <- sprintf("F=%s", config.forecast.scenario1.f_year1)
name(result.forecast.F_scenario2) <- sprintf("F=%s", config.forecast.scenario2.f_year1)
name(result.forecast.F_SQ) <- sprintf("F(SQ)=%s", round(F_SQ, 2))

# compile summary object with all forecasts
result.forecast.all <- FLStocks(result.forecast.F_scenario1, result.forecast.F_scenario2, result.forecast.F_SQ, result.forecast.F_01manual, result.forecast.F_MSY)

# save output
savePic("output/8_Forecasts/summary_all_in_one.png", result.forecast.all)
savePic("output/8_Forecasts/short_all_in_one.png", window(result.forecast.all, start=lastYear))


print("[XSA toolkit]: Procedure Done! Take a look on /output/ folder")

render("Result.Rmd")

